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Abstract The resistive or non-resistive nature of the 
extracellular space in the brain is still debated, and is an 
important issue for correctly modeling extracellular 
potentials. Here, we first show theoretically that if the 
medium is resistive, the frequency scaling should be the 
same for electroencephalogram (EEG) and 
magnetoencephalogram (MEG) signals at low frequencies 
(<10 Hz). To test this prediction, we analyzed the 
spectrum of simultaneous EEG and MEG measurements 
in four human subjects. The frequency scaling of EEG 
displays coherent variations across the brain, in general 
between 1/f and l/f 2 , and tends to be smaller in 
parietal/temporal regions. In a given region, although the 
variability of the frequency scaling exponent was higher 
for MEG compared to EEG, both signals consistently 
scale with a different exponent. In some cases, the scaling 
was similar, but only when the signal-to-noise ratio of the 
MEG was low. Several methods of noise correction for 
environmental and instrumental noise were tested, and 
they all increased the difference between EEG and MEG 
scaling. In conclusion, there is a significant difference in 
frequency scaling between EEG and MEG, which can be 
explained if the extracellular medium (including other 
layers such as dura matter and skull) is globally 
non-resistive. 
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1 Introduction 



An issue central to modeling local field potentials is 
whether the extracellular space around neurons can be 
considered as a resistive medium. A resistive medium is 
equivalent to replacing the medium by a simple 
resistance, which considerably simplifies the computation 
of local field potentials, as the equations to calculate 
extracellular fields are very simple and based on 
Coulomb's law (Rail and Shepherd, 1968; Nunez and 
Srinivasan, 2005). Forward models of the EEG and 
inverse solution/source localization methods also assume 
that the medium is resistive (Sarvas, 1987; Wolters and de 
Munck, 2007; Ramirez, 2008). However, if the medium is 
non-resistive, the equations governing the extracellular 
potential can be considerably more complex because the 
quasi-static approximation of Maxwell equations cannot 
be made (Bedard et al., 2004). 

Experimental characterizations of extracellular resistivity 
are contradictory. Some experiments reported that the 
conductivity is strongly frequency dependent, and thus 
that the medium is non-resistive (Ranck, 1963; Gabriel et 
al., 1996a, 1996b, 1996c). Other experiments reported 
that the medium was essentially resistive (Logothetis et 
al., 2007). However, both types of measurements used 
current intensities far larger than physiological currents, 
which can mask the filtering properties of the tissue by 
preventing phenomena such as ionic diffusion (Bedard 
and Destexhe, 2009). Unfortunately, the issue is still open 
because there exists no measurements to date using 
(weak) current intensities that would be more compatible 
with biological current sources. 

In the present paper, we propose an indirect method to 
estimate if extracellular space can be considered as a 
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purely resistive medium. We start from Maxwell 
equations and show that if the medium was resistive, the 
frequency-scaling of electroencephalogram (EEG) and 
magnetoencephalogram (MEG) recordings should be the 
same. We then test this scaling on simultaneous EEG and 
MEG measurements in humans. 



2 Methods 

2.1 Participants and MEG/EEG recordings 

We recorded the electromagnetic field of the brain during 
quiet wakefulness (with alpha rhythm occasionally 
present) from four healthy adults (4 males ages 20-35). 
Participants had no neurological problems including sleep 
disorders, epilepsy, or substance dependence, were taking 
no medications and did not consume caffeine or alcohol 
on the day of the recording. We used a whole-head MEG 
scanner (Neuromag Elekta) within a magnetically 
shielded room (IMEDCO, Hagendorf, Switzerland) and 
recorded simultaneously with 60 channels of EEG and 
306 MEG channels (Nenonen et al., 2004). MEG SQUID 
(super conducting quantum interference device) sensors 
are arranged as triplets at 102 locations; each location 
contains one "magnetometer" and two orthogonal planar 
"gradiometers" (GRAD1, GRAD2). Unless otherwise 
noted, MEG will be used here to refer to the 
magnetometer recordings. Locations of the EEG 
electrodes on the scalp of individual subjects were 
recorded using a 3D digitizer (Polhemus FastTrack). HPI 
(head position index) coils were used to measure the 
spatial relationship between the head and scanner. 
Electrode arrangements were constructed from the 
projection of 3D position of electrodes to a 2D plane in 
order to map the frequency scaling exponent in a 
topographical manner. All EEG recordings were 
monopolar with a common reference. Sampling rate was 
1000 Hz. 

For all subjects, four types of consecutive recordings were 
obtained, in the following order: (1) Empty -room 
recording; (2) Awake "idle" recording where subjects 
were asked to stay comfortable, without movements in the 
scanner, and not to focus on anything specific; (3) a visual 
task; (4) sleep recordings. All idle recordings used here 
were made in awake subjects with eyes open, where the 
EEG was desynchronized. A few minutes of such idle 
time was recorded in the scanner. For each subject, 3 
awake segments with duration of 60 seconds were 
selected from the idle recordings (see example signals in 
Fig.ID. 



As electrocardiogram (ECG) noise often contaminates 
MEG recordings, Independent component analysis (ICA) 
algorithm was used to remove such contamination; either 
Infomax (Bell and Sejnowski, 1995) or the "Jade 
algorithm" from the EEGLAB toolbox (Delorme and 
Makeig, 2004) was used to achieve proper 
decontamination. In all recordings, the ECG component 
stood out very robustly. In order not to impose any change 
in the frequency content of the signal, we did not use the 
ICA to filter the data on any prominent independent 
oscillatory component and it was solely used to 
decontaminate the ECG noise. We verified that the 
removal of ECG did not change the scaling exponent (not 
shown). 

In each recording session, just prior to brain recordings, 
we recorded a few minutes of the electromagnetic field 
present within the dewar in the magnetic shielded room. 
Similar to wake epochs, 3 segments of 60 seconds 
duration were selected for each of the four recordings. 
This will be referred to "empty room" recordings and will 
be used in noise correction of the awake recordings. 

In each subject, the power spectral density (PSD) was 
calculated by first computing the Fast Fourier transform 
(FFT) of 3 awake epochs, then averaging their respective 
PSDs (square modulus of the FFT). This averaged PSD 
was computed for all EEG and MEG channels in order to 
reduce the effects of spurious peaks due to random 
fluctuations. The same procedure was also followed for 
empty-room signals. 



2.2 Noise correction methods 

Because the environmental and instrumental sources of 
noise are potentially high in MEG recordings, we took 
advantage of the availability of empty-room recordings to 
correct for the presence of noise in the signal. We used 
five different methods for noise correction, based on 
different assumptions about the nature of the noise. We 
describe below these different correction methods, while 
all the details are given in Supplementary Methods. 

A first procedure for noise correction, exponent 
subtraction (ES), assumes that the noise is intrinsic to the 
SQUID sensors. This is justified by the fact that the 
frequency scaling of some of the channels is identical to 
that of the corresponding empty-room recording (see 
Results). In such a case, the scaling is assumed to entirely 
result from the "filtering" of the sensor, and thus the 
correction amounts to subtract the scaling exponents. 



3 



A second class of noise subtraction methods assume that 
the noise is of ambient nature and is uncorrected with the 
signal. This chatacteristics, warrants the use of spectral 
subtraction (where one subtracts the PSD of the 
empty-room from that of the MEG recordings), prior to 
the calculation of the scaling exponent. The simplest form 
of spectral subtraction, linear multiband spectral 
subtraction (LMSS), treats the sensors individually and 
does not use any spatial/frequency-based statistics in its 
methodology (Boll et al., 1979). An improved version, 
nonlinear multiband spectral subtraction (NMSS), takes 
into account the signal-to-noise ratio (SNR) and its spatial 
and frequency characteristics (Kamath and Loizou, 2002; 
Loizou, 2007). A third type, Wiener filtering (WF), uses a 
similar approach as the latter, but obtain an estimate of the 
noiseless signal from that of the noisy measurement 
through minimizing the Mean Square Error (MSE) 
between the desired and the measured signal (Lim et al., 
1979; Abd El-Fattah et al., 2008). 

A third type of noise subtraction, partial least squares 
(PLS) regression, combines Principal component analysis 
(PCA) methods with multiple linear regression (Abdi, 
2010; Garthwaite, 1994). This methods finds the spectral 
patterns that are common in the MEG and the 
empty-room noise, and removes these patterns from the 
PSD. 

2.3 Frequency scaling exponent estimation 

The method to estimate the frequency scaling exponent 
was composed of steps: First, applying a spline to obtain a 
smooth FFT without losing the resolution (as can happen 
by using other spectral estimation methods); Second, 
using a simple polynomial fit to obtain the scaling 
exponent. To improve the slope estimation, we 
approximated the PSD data points using a spline, which is 
a series of piecewise polynomials with smooth transitions 
and where the break points ("knots") are specified. We 
used the so-called "B-spline" (see details in de Boor, 
2001). 

The knots were first defined as linearly related to 
logarithm of the frequency, which naturally gives more 
resolution to low frequencies, to which our theory applies. 
Next, in each frequency window (between consecutive 
knots), we find the closest PSD value to the mean PSD of 
that window. Then we use the corresponding frequency as 
the optimized knot in that frequency range, leading the 
final values of the knots. The resulting knots stay close to 
the initial distribution of frequency knots but are modified 
based on each sensor's PSD data to provide the optimal 
knot points for that given sensor (Fig. EJV). We also use 
additional knots at the outer edges of the signal to avoid 



boundary effects (Eilers and Marx, 1996). The applied 
method provides a reliable and automated approach that 
uses our enforced initial frequency segments with a high 
emphasis in low frequency and it optimizes itself based on 
the data. After obtaining a smooth B-spline curve, a 
simple 1st degree polynomial fit was used to estimate the 
slope of the curve between 0.1-10 Hz (the fit was limited 
to this frequency band in order to avoid the possible 
effects of the visible peak at 10 Hz on the estimated 
exponent). Using this method provides a reliable and 
robust estimate of the slope of the PSD in logarithmic 
scale, as shown in Fig. |2]B. For more details on the issue 
of automatic non-parametric fitting, and the rationale 
behind combining the polynomial with spline basis 
functions, we refer the reader to Magee, 1998 as well as 
Royston & Altman, 1994 and Katkovnik et al, 2006. 

This procedure was realized on all channels automatically 
(102 channels for MEG, 60 channels for EEG, for each 
patient). Every single fit was further visually confirmed. 
In the case of MEG, noise correction is essential to 
validate the results. For doing so, we used different 
methods (as described above) to reduce the noise. Next, 
all the mentioned steps of frequency scaling exponents 
were carried out on the corrected PSD. Results are shown 
in Fig. E] 



2.4 Region of Interest (ROI) 

Three ROIs were selected for statistical comparisons of 
the topographic plots. As shown in Figure [4] (panel F), FR 
(Frontal) ROI refers to the frontal ellipsoid, VX (Vertex) 
ROI refers to the central disk located on vertex and PT 
(Parietotemporal) refers to the horseshoe ROI. 



3 Theory 

We start from first principles (Maxwell equations) and 
derive equations to describe EEG and MEG signals. Note 
that the formalism we present here is different than the 
one usually given (as in Plonsey, 1969; Gulrajani, 1998), 
because the linking equations are here considered in their 
most general expression (convolution integrals), in the 
case of a linear medium (see Eq. 77.4 in Landau and 
Lifchitz, 1984). This generality is essential for the 
problem we treat here, because our aim is to compare 
EEG and MEG signals with the predictions from the 
theory, and thus the theory must be as general as possible. 
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3.1 General formalism 



Maxwell equations can be written as 



VD = p 



Jree 



VxE=-^ 



V-B = 
VxH = j + f» 



(1) 



By taking the Fourier transform of Maxwell equations 
(Eqs. [TJ and of the linking equations (Eqs. 121314b . we 
obtain: 
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(6) 



If we suppose that the brain is linear in the 
electromagnetic sense (which is most likely), then we 
have the two following linking equations. The first 
equation links the electric displacement with the electric 
field: 



D 



e{i)E{t ~ x)dx 



(2) 



where e is a symmetric second-order tensor. 



A second equation links magnetic induction and the 
magnetic field: 



B 



(3) 



where jj, is a symmetric second-order tensor. 



If we neglect non-resistive effects such as diffusion 
(Bedard and Destexhe, 2009), as well as any other 
nonlinear effectfl then we can assume that the medium is 
linear. In this case, we can write: 



/ — ^/ / 
if =f f + o f E f 



(7) 



where the relation oyE t in Eq. [7] is the current density 
produced by the (primary) current sources in the 
extracellular medium. Note that in this formulation, the 
electromagnetic parameters £/, jif and oy depend on 
frequency This generalization is essential if we want the 
formalism to be valid for media that are linear but 
non-resistive, which can expressed with 
frequency-dependent electric parameters. It is also 
consistent with the Kramers-Kronig relations (see Landau 
and Lifchitz, 1984; Foster and Schwan, 1989). 

is the current density of these sources in Fourier 
frequency space. This current density is composed of the 
axial current in dendrites and axons, as well as the 
transmembrane current. Of course, this expression is such 
that at any given point, there is only one of these two 
terms which is non-zero. This is a way of preserving the 
linearity of Maxwell equations. Such a procedure is 
legitimate because the sources are not affected by the field 
they product 



J 



C7(T)E(f-T)rfT 



(4) 



where a is a symmetric second-order tensor^. Because the 
effect of electric induction (Faraday's law) is negligible, 
we can write: 



V D = p free 



V B = 



VxE 



VxH = j 



<3D 

~3F 



(5) 



This system is much simpler compared to above, because 
electric field and magnetic induction are decoupled. 

1 Examples of nonlinear effects are variations of the macroscopic 
conductivity 07 with the magnitude of electric field E. Such 
variations could appeal' due to ephaptic (electric-field) interactions 
for example. In addition, any type of linear reactivity of the medium 
to the electric field or magnetic induction can lead to frequency- 
dependent electric parameters <T,e,/i (for a detailed discussion of 
such effects, see Bedard and Destexhe, 2009). 

2 Note that in textbooks, these linking equations (Eqs. are 
often algebraic and independent of time (for example, see Eqs. 5.2-6, 
5.2-7 and 5.2-8 in Gulrajani, 1998). The present formulation is more 
general, more in the line of Landau and Lifchitz (1984). 



3.2 Expression for the electric field 



From Eq. [6](Faraday's law in Fourier space), we can write: 



E f = -W, 



(8) 



From Eq.[6]( Ampere-Maxwell's law in Fourier space), we 
can write: 



V • (V x H/) = V-j f + icoV ■ (SfEf) 



V -j/-V -((Of + i(oe f ) VV f ) =0 



(9) 



3 In textbooks, the electric parameters are sometimes considered 
as complex numbers, for example with the notion of phasor (see 
Section 5.3 in Gulrajani, 1998), but they are usually considered 
frequency independent. 

4 If it was not the case, then the source terms would be a function 
of the produced field, which would result in more complicated 
equations 
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Setting Yf = Of + icaEf, one obtains: 
V-(Y f VVf)=V-f f 



(10) 



where V ■ is a source term and yy is a symmetric 
second-order tensor (3 x 3). Note that this tensor depends 
on position and frequency in general, and cannot be 
factorized. We will call this expression (Eq. ITOb the "first 
fundamental equation" of the problem. 



3.3 Expression for magnetic induction 



From the mathematical identity 
Vx VxX=-V 2 X + V(V-X) 



(11) 



it is clear that this is sufficient to know the divergence and 
the curl of a field X, because the solution of V 2 X is unique 
with adequate boundary conditions. 

As in the case of magnetic induction, the divergence is 
necessarily zero, it is sufficient to give an explicit 
expression of the curl as a function of the sources. 

Supposing that ji = ji 8(t) is a scalar (tensor where all 
directions are eigenvectors), and taking the curl of Eq. [6] 
(D), multiplied by the inverse of Yf> we obtain the 
following equality: 



(12) 



V x (yf'V x B f ) = H V x (y}%) 

because V x Ey = 0. This expression (Eq. [T2l will be 
named the "second fundamental equation". 



3.4 Boundary conditions 



We consider the following boundary conditions: 

1 - on the skull, we assume that Vy(r) is differentiable in 
space, which is equivalent to assume that the electric field 
is finite. 

2 - on the skull, we assume that n ■ yyVVy is also 
continuous, which is equivalent to assume that the flow of 
current is continuous. Thus, we are interested in solutions 
where the electric field is continuous. 

3 - because the current is zero outside of the head, the 
current perpendicular to the surface of cortex must be zero 
as well. Thus, the projection of the current on the vector h 
normal to the skull's surface, must also be zero. 



The latter expression can be proven by calculating the total 
current and apply the divergence theorem (not shown). 



3.5 Quasi-static approximation to calculate magnetic 
induction 



The "second fundamental equation" above implies 
inverting Yf> which is not possible in general, because it 
would require prior knowledge of both conductivity and 
permittivity in each point outside of the sources. If the 
medium is purely resistive (yy = y where y is independent 
of space and frequency), one can evaluate the electric field 
first, and next integrate By using the quasi-static 
approximation (Ampere-Maxwell's law). Because for low 
frequencies, we have necessarily j f >> /a)Dy, we obtain 

VxB / = ju j / , 

which is also known as Ampere's law in Fourier space. 

Thus, for low frequencies, one can skip the second 
fundamental equation. Note that in case this quasi-static 
approximation cannot be made (such as for high 
frequencies), then one needs to solve the full system using 
both fundamental equations. Such high frequencies are, 
however, well beyond the physiological range, so for EEG 
and MEG signals, the quasi-static approximation holds if 
the extracellular medium is resistive, or more generally if 
the medium satisfies V x Ey = — /©By ^ (see Eqs. [5] 
and|6]l. 

According to the quasi-static approximation, and using the 
linking equation between current density and the electric 
field (Eq.|7]i, we can write: 



(14) 



Because the divergence of magnetic induction is zero, we 
have from Eq. QT| 



V x V x By - -V 2 By = -jfcV x (j'; - yWy) 



(15) 



This equation can be easily integrated using Poisson 
integral ("Poisson equation" for each component in 
Cartesian coordinates) In Fourier space, this integral is 
given by the following expression 



Vx(jP(r')-yWy(r')) 



dv' 



(16) 



head 



n(x) ■ yyVVy(x) = 



(13) 
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3.6 Consequences 



If the medium is purely resistive ("ohmic"), then y does 
not depend on the spatial position (see Bedard et al., 2004; 
Bedard and Destexhe, 2009) nor on frequency, so that the 
solution for the magnetic induction is given by: 

head 

and does not depend on the nature of the medium. 

For the electric potential, from Eq. [10] we obtain the 
solution: 



V f (x) 



47ZJ JJJ |x-x' 

head 



dv' 



(18) 



Thus, when the two source terms V x and V • are 
white noise, the magnetic induction and electric field must 
have the same frequency dependence. Moreover, because 
the spatial dimensions of the sources are very small (see 
appendices), we can suppose that the current density j^(x) 
is given by a function of the form: 



j?(x)=r(x)F(/) 



(19) 



such that V x and V ■ have the same frequency 
dependence for low frequencies. Eq. [19] constitutes the 
main assumption of this formalism. 

In Appendix A, we provide a more detailed justification 
of this assumption, based on the differential expressions 
of the electric field and magnetic induction in a dendritic 
cable. Note that this assumption is most likely valid for 
states with low correlation such as desynchronized-EEG 
states or high-conductance states, and for low-frequencies, 
as we analyze here (see details in the appendices). 

Thus, the main prediction of this formalism is that if the 
extracellular medium is resistive, then the PSD of the 
magnetic induction and of the electric potential must have 
the same frequency dependence. In the next section, we 
will examine if this is the case for simultaneously 
recorded MEG and EEG signals. 



had eyes open, a low-amplitude alpha rhythm was 
occasionally present (as visible in Fig. [TJ. There were also 
oscillations present in the empty-room signal, but these 
oscillations are evidently different from the alpha rhythm 
because of their low amplitude and the fact that they do 
not appear in gradiometers (see Suppl. Fig. SI). 



A EEG 



B 



MEG 



1 sec 



******* 



Fig. 1 Simultaneous EEG and MEG recordings in an awake human 
subject. This example shows a sample of channels from MEG/EEG 
after ECG noise removal. Labels refer to ROIs as defined in methods 
(also see Figure|4j. FR: Frontal, VX:Vertex and PT: Parietotemporal. 
These sample channels were selected to represent both right and left 
hemispheres in a symmetrical fashion. Inset: magnification of the 
MEG (red) and "empty-room" (green) signals superimposed from 4 
sample channels. All traces are before any noise correction, but after 
ECG decontamination. 



In the next sections, we start by briefly presenting the 
method that was used to estimate the frequency scaling of 
the PSDs. Then we report the scaling exponents for 
0.1-10 Hz frequency bands and their differences in EEG 
and MEG recordings. 



4. 1 Frequency scaling exponent estimation 



4 Test on experimental data 

A total of 4 subjects were used for the analysis. Figure Q] 
shows sample MEG and EEG channels from one of the 
subjects, during quiet wakefulness. Although the subjects 



Because of the large number of signals in the EEG and 
MEG recordings, we used an automatic non-parametric 
procedure to estimate the frequency scaling (see 
Methods). We used a B-spline approximation by 
interpolation with boundary conditions to find a curve 
which best represents the data(see Methods). A high 
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density of knots was given to the low-frequency band 
(0.1-10 Hz), to have an accurate representation of the PSD 
in this band, and calculate the frequency scaling. An 
example of optimized knots to an individual sensor is 
shown in Figure [2|\; note that this distribution of knots is 
specific to this particular sensor. The resulting B-spline 
curves were used to estimate the frequency scaling 
exponent using a 1st degree polynomial fit. Figure |2jj 
shows the result of the B-spline analysis with optimized 
knots (in green) capturing the essence of the data better 
than the usual approximation of the slope using 
polynomials (in red). The goodness of fit showed a robust 
estimation of the slope using B-spline method. Residuals 
were -0.01 ± 0.6 for empty-room, 0.2 ± 0.65 for MEG 
awake, 0.05 ± 0.6 for LMSS, 0.005 ± 0.64 for NMSS, 
0.08 ± 0.5 for WF,0.001 ± 0.02 for PLS, and -0.02 ± 
0.28 for EEG B-spline (all numbers to be multiplied by 

io- 14 ). 
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Fig. 2 A.log-log scale of the PSD vs frequency of a sample MEG 
sensor along with the corresponding log(PSD) values (shown as 
circles) at optimized knots in log-scale. B. 1st degree Polynomial fit 
on B-spline curve effectively captures properties of the signal better 
than simple polynomial fit and avoids the 10 Hz peak. The fit was 
limited between 0.1 to 10 Hz excluding the boundaries. This limits 
the fit approximation to the next limiting optimized knots (between 
0.1 and 0.2 to between 9 and 10 Hz) to avoid the peaks at alpha and 
low frquencies (shown by vertical dotted lines). 





— B-spline of PSD 




— 1st orderjpoly fit to PSD 




>i^— 1st orderjpoly fit to B-spline 






2 -1 


1 2 




Log(f) 



4.2 MEG and EEG have different frequency scaling 
exponents 

Figure [3] shows the results of the B-spline curve fits to the 
log-log PSD vs frequency for all sensors of all subjects. 
In this figure, and only for the ease of visual comparison, 
these curves were normalized to the value of the log(PSD) 
of the highest frequency. As can be appreciated, all MEG 
sensors (in red) show a different slope than that of the EEG 
sensors (in blue). The frequency scaling exponent of the 
EEG is close to 1 (1 // scaling), while MEG seems to scale 
differently. Thus, this representation already shows clear 
differences of scaling between EEG and MEG signals. 



Subject 1 



Subject 2 




Subject 3 




Subject 4 




Fig. 3 B-spline fits of EEG awake and MEG awake (prior to noise 
correction) recordings from all four subjects. Each line refers to the fit 
of one sensor in log(PSD)-log(frequency) scale. For the ease of visual 
comparison of the frequency scaling exponent, log(PSD) values are 
normalized to their value at the maximum frequency. Each panel 
represents the data related to one of our four subjects. These plots 
show a clear distinction between the frequency scaling of EEG and 
MEG. Insets show the comparison between MEG awake (prior to 
noise correction) and MEG empty-room recordings (not normalized). 
Note that the empty-room scales the same as the MEG signal, but in 
general EEG and MEG scale differently. 



However, MEG signals may be affected by ambient or 
instrumental noise. To check for this, we have analyzed 
the empty-room signals using the same representation and 
techniques as for MEG, amd the results are represented in 
Fig. [3] (insets). Empty-room recordings always scale very 
closely to the MEG signal, and thus the scaling observed 
in MEG may be due in part to environmental noise or 
noise intrinsic to the detectors. This emphasizes that it is 
essential to use empty-room recordings made during the 
same experiment to correct the frequency scaling 
exponent of MEG recordings. 

To correct for this bias, we have used five different 
procedures (see Methods). The first class of procedure 
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(ES) considers that the scaling of the MEG is entirely due 
to filtering by the sensors, which would explain the 
similar scaling between MEG and empty-room 
recordings. In this case, however, nearly all the scaling 
would be abolished, and the corrected MEG signal would 
be similar to white noise (scaling exponent close to zero). 
Because the similar scaling may be coincidental, we have 
used two other classes of noise correction procedures to 
comply with different assumptions about the nature of the 
noise. The second class, is composed of spectral 
subtraction (LMSS and NMSS) or Wiener filtering (see 
Methods). These methods are well-established in other 
fields such as acoustics. The third class, uses statistical 
patterns of noise to enhance PSD (PLS method, for details 
see Methods). 



4.3 Spatial variability of the frequency scaling exponent 

We applied the above methods to all channels and 
represented the scaling exponents in topographic plots in 
Fig. |4] This figure portrays that both MEG and EEG do 
not show a homogenous pattern of the scaling exponent, 
confirming the differences of scaling seen in Fig. [3] The 
EEG (Figure |4j\) shows that areas in the midline have 
values closer to 1, while those at the margin can deviate 
from l/f scaling. MEG on the other hand shows higher 
values of the exponent in the frontal area and a horseshoe 
pattern of low value exponents in parietotemporal regions 
(Figure [4)3). As anticipated above, empty-room 
recordings scale more or less uniformly with values close 
to l/f (Figure Hf), thus necessitating the correction for 
this phenomena to estimate the correct MEG frequency 
scaling exponent. Different methods for noise reduction 
are shown in Figure |4j spectral subtraction methods, such 
as LMSS (Figure gt>), NMSS (Figure HE), WF 
enhancement (Figure Hf). These corrections preserve the 
pattern seen in Figure |4]B, but tend to increase the 
difference with EEG scaling: one method (LMSS) yields 
minimal correction while the other two (NMSS and WF) 
use band-specific SNR information in order to cancel the 
effects of background colored-noise (see Suppl. Fig. S2), 
and achieve higher degree of correction (see 
Supplementary Methods for details). Figure HJj portrays 
the use of PLS to obtain a noiseless signal based on the 
noise measurements. The degree of correction achieved 
by this method is higher than what is achieved by spectral 
subtraction and WF methods. Exponent subtraction is 
shown in Figure |4jL This correction supposes that the 
scaling is due to the frequency response of the sensors, 
and nearly abolishes all the frequency scaling (see also 
Suppl. Fig. S3 for a comparison of different methods of 
noise subtraction). 




Fig. 4 Topographical representation of frequency scaling exponent 
averaged across four subjects. A. EEG awake. B. MEG awake. C. 
MEG empty-room. D, E. MEG after spectral subtraction of the 
empty-room noise using linear (LMSS) and non-linear (NMSS) 
methods respectively. F. MEG spectral enhancement using Wiener 
filtering (WF). G. MEG, partial least square (PLS) approximation 
of non-noisy spectrum. H. Exponent subtraction (the exponent 
represented is the value of the frequency scaling exponent calculated 
for MEG signals, subtracted from the scaling exponent calculated 
from the corresponding emptyroom signals). I. Spatial location of 
ROI masks (shown in yellow). FR covers the Frontal, VX covers 
Vertex and PT spans Parietotemporal. Dots show spatial arrangement 
of 102 MEG SQUID sensor triplets. The background gray-scale 
figure is same as the one in panel B. Note that panels A through H 
use the same color scaling. 



4.4 Statistical comparison of EEG and MEG frequency 
scaling 

Based on the patterns in Fig. [4] we created three ROIs 
covering Vertex (FR), Vertex (VX) and the horseshoe 
pattern (PT). These masks are shown in Fig.|4j. 

Figure |5}\ represents the overall pattern providing 
evidence on the general difference and the wider 
variability in MEG recordings. The next three panels 
relate to the individual ROIs. Of the spectral subtraction 
methods, NMSS achieves a higher degree of correction in 
comparison with LMSS (see Figure[4f , Figure|4]D as well 
as Suppl. Fig. S3). Because NMSS takes into account the 
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A. All areas 




-2 -1.5 -1 -0.5 0.5 

D. Parietotemporal (PT) ROI 



NMSS 
MEG 




0.05 





Fig. 5 Statistical comparison of EEG vs. MEG frequency scaling 
exponent for all regions (A) and different ROI masks (B,C & 
D). In each panel, a box-plot on top is accompanied by a non- 
parametric distribution function in the bottom. In the top graph, the 
box has lines at the lower quartile, median (red), and upper quartile 
values. Smallest and biggest non-outlier observations (1.5 times the 
interquartile range IRQ) are shown as whiskers. Outliers are data with 
values beyond the ends of the whiskers and are displayed with a 
red + sign. In the bottom graph, a Non-parametric density function 
shows the distribution of EEG, MEG and empty-room-corrected 
MEG frequency scaling exponents (note that LMSS and WF are not 
shown here; see the text for description.). Thick and thin vertical lines 
show the mean and mean ± std for each probability density function 
(pdf). 



effects of the background colored-noise (Suppl. Fig. S2), 
it is certainly more relevant to the type of signals analyzed 
here. The results of NMSS and WF are almost identical 
and confirm one another (see Figure |4js, as well as Suppl. 



Fig. S3). Therefore, of this family of noise correction, 
only NMSS is portrayed here. Of the methods dealing 
with different assumptions about the nature of the noise, 
the "Exponent subtraction" almost abolishes the 
frequency scaling (Also see in Figure |4}L as well as 
Suppl. Fig. S3). Applying PLS yields values in between 
"Exponent subtraction" and that of NMSS and is 
portrayed in Figure [5] 

In the Frontal region (Figure |5j3), the EEG scaling 
exponents show higher variance by comparison to MEG. 
Also, EEG shows some overlaps with the distribution 
curve of non-corrected MEG; this overlap becomes 
limited to the tail end of the NMSS correction and is 
abolished in the case of PLS correction. As can be 
appreciated, VX (Figure |5p) shows both similar values 
and similar distribution for EEG and non-corrected MEG. 
These similarities, in terms of regional overall values and 
distribution curve, are further enhanced after NMSS 
correction. It is to be noted that, in contrast to these 
similarities, the one-to-one correlation of NMSS and EEG 
at VX ROI are very low (see below, Table 1B-C). The 
values of PLS noise correction are very different from that 
of EEG and have a similar, but narrower, distribution 
curve shape. Two other ROIs show distinctively different 
values and distribution in comparing EEG and MEG. 
Both NMSS and PLS agree on this with PLS showing 
more extreme cases. Figure |5f) reveals a bimodal 
distribution of MEG exponents in the parietotemporal 
region (PT ROI). This region has also the highest variance 
(in MEG scaling exponents) compared to other ROIS. The 
distinction between EEG and MEG is enhanced in PLS 
estimates; however, the variance of PT is reduced in 
comparison to NMSS while the bimodality is still 
preserved but weakened. The values of mean and standard 
deviation for these ROIs' exponents are provided in 
Table 1 A (mean ± standard deviation). 

The box-plots of Fig. [5}plots further show the difference 
between the medians, lower/upper quartile and 
interquartile range. The overall difference is that the 
uncorrected MEG has much wider variance compared to 
EEG and corrected MEG (in case of PLS correction); the 
absolute value of the median of MEG (uncorrected, or 
corrected with either NMSS or PLS) is always smaller 
than that of EEG. The VX region is an exception to the 
above rules; interestingly, the one-to-one correlation of 
VX happens to be the lowest of all (see below). In the 
case of NMSS -corrected MEG, the shape of the pdf is 
preserved. However, PLS narrows the distribution curve 
of MEG but further enhances the differences between 
MEG and EEG. Therefore, median and lower/upper 
quartiles will have different value than that of EEG. 
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A. Mean and standard deviation 





EEC 




All 


-1.33 ±0.19 


- 1 24 ± 26 - 1 06 ± 29 


FRROI 


-1.36 ±0.25 


-0.97 ±0.10 -0.76 ±0.09 


VXROI 


-1.21 ±0.13 


-1.36 ±0.10 -1.14 ±0.11 


PTROI 


-1.36 ±0.12 


-1.30 ±0.29 -1.16 ±0.32 


B. Pearson correlation 






EEG vs. MEG 


EiCVJ VS. v-OI 1 CL LcU 1V1CU ^INlVlOO^ 


All 


0.29 


0.32 


FRROI 


0.41 


0.32 


VXROI 


-0.17 


-0.15 


PTROI 


0.35 


0.38 


C. Kendall Rank Corr 






EEG vs. MEG 


EEG vs Corrected MEG (NMSS) 


All 


0.21 


0.24 


FRROI 


0.29 


0.21 


VXROI 


-0.03 


-0.04 


PTROI 


0.23 


0.26 



D. KruskalWallis 





p value 


Chi-square 


df Error 


All 


< io- b 


1.53 10 3 


34838 


All noise-corrected 


< io- 15 


8.03 10 3 


34838 


FRROI 


< io- 15 


3.30 10 3 


5008 


FR ROI noise-corrected 


< io- 15 


3.72 10 3 


5008 


VXROI 


< io- 15 


1.72 10 3 


5452 


VX ROI noise-corrected 


< io- 15 


0.23 10 3 


5452 


PTROI 


< io- 15 


0.21 10 3 


13010 


PT ROI noise-corrected 


< io- 15 


1.18 10 3 


13010 


Table 1 ROI statistical 


comparison. 


A. mean 


and std 



frequency scale exponent for all regions and individual ROI. 
B. numerical values of linear Pearson correlation. C. rank-based 
Kendall correlation. D. non-parametric test of analysis of variance 
(KruskalWallis). Corrected MEG refers to spectral subtraction using 
NMSS. The full table is provided in Supplementary information. 



Correlation values (Table 1B-C) show that, although VX 
ROI has the closest similarity in terms of its central 
tendency and probability distribution, it provides the 
lowest correlation in a pairwise fashion. P-values (for 
testing the hypothesis of no correlation against the 
alternative that there is a nonzero correlation) for 
Pearson's correlation were calculated using a Student's 
t-distribution for a transformation of the correlation and 
they were all significant (less than 10~ 15 for a = 0.05). 
Similarly, a non-parametric statistic Kendall tau rank 
correlation was used to measure the degree of 
correspondence between two rankings and assessing the 
significance of this correspondence between MEG and 
EEG in the selected ROIs (Table 1C). P-values for 
Kendall's tau and Spearman's rho calculate using the 
exact permutation distributions were all significant (less 
than 10~ 15 for a = 0.05). Kendall tau shows that the rank 
correlation for all areas considered together as well as for 
PT, show a lesser correlation than that is shown by 
Pearson linear correlation. Furthermore, we carried out a 
Kruskal-Wallis nonparametric version of one-way 
analysis of variance. We used this test to avoid bias in 
ANOVA (KruskalWallis assumes that the measurements 
come from a continuous distribution, but not necessarily a 



normal distribution as is assumed in ANOVA). 
KruskalWallis uses analysis of variance on the ranks of 
the data values, not the data values themselves and 
therefore is an appropriate test for comparison of the 
homogeneity of pattern between ROIs of two image as 
well as their statistical median. As shown in Table ID, all 
p-values were significant emphasizing the difference 
between the spatial aspect of the spectral nature of MEG 
and EEG. Note that the difference of scaling exponent of 
EEG and MEG was also confirmed by nonlinear spatial 
kendall correlation analysis, independently of the ROIs 
classification (not shown). 



4.5 Relation of scaling exponent to signal-to-noise ratio 

Noise correction does not affect all the sensors in a same 
fashion. As presented in Suppl. Fig. S3, the simple linear 
spectral subtraction (LMSS) may lead to an increment or 
decrement of the scaling exponent. In any case, the 
correction achieved by this method is minimal. This is due 
to the fact that LMSS ignores the complex non-linear 
patterns of the SNR in different channels (Suppl. Fig. S2). 
We show that for all subjects, as the frequency goes up, 
the SNR goes down. It is also noticable that in each 
defined frequency band, i.e. 0-10 Hz (Slow, Delta and 
Theta), 11-30 Hz (Beta), 30-80 Hz (Gamma), 80-200 Hz 
(Fast oscillation), 200-500 Hz (Ultra-fast oscillation), 
there is an observable sensor-to-sensor SNR variability. 
This variability is at its maximum in the band with the 
highest SNR (i.e. 1-10 Hz). All together, the non-linear 
nature of MEG SNR shows that a linear spectral 
subtraction could behave non-optimally, leading to 
minimal correction. This also conveys that the optimal 
spectral correction can be achieved only by non-linear 
methods that explicitly take into account the SNR 
information of the data. Therefore the correction achieved 
by NMSS and WF have higher validity, in agreement with 
the fact that both methods yield similar results in terms of 
values and spatial distribution (Fig. 2)2, Fig. Hp). 

5 Discussion 

In this paper, we have used a combination of theoretical 
and experimental analyses to investigate the spectral 
structure of EEG and MEG signals. In the first part of the 
paper, we presented a theoretical investigation showing 
that if the extracellular medium is purely resistive, the 
equations of the frequency dependence of electric field 
and magnetic induction take a simple form, because the 
admittance tensor does not depend on spatial coordinates. 
Thus, the macroscopic magnetic induction does not 
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depend on the electric field outside the neuronal sources, 
but only depends on currents inside neurons. In this case, 
the frequency scaling of the PSD should be the same for 
EEG and MEG signals. This conclusion is only valid in 
the linear regime, and for low frequencies. 

An assumption behind this formalism is that the spatial 
and frequency dependence of the current density factorize 
(Eq. |T~9T >. We have shown in the appendices that this is 
equivalent to consider the different current sources as 
independent. Thus, the formalism will best apply to states 
where the activity of synapses is intense and of very low 
correlation. This is the case for desynchronized-EEG 
states or more generally "high-conductance states", in 
which the activity of neurons is intense, of low 
correlation, and the neuronal membrane is dominated by 
synaptic conductances (Destexhe et al., 2003). In such 
conditions, the dendrites are bombarded by intense 
synaptic inputs which are essentially uncorrelated, and 
one can consider the current sources as independent 
(Bedard et al., 2010). In the present paper, we analyzed 
EEG and MEG recordings in such desynchronized states, 
where this formalism best applies. 

Note that the above reasoning neglects the possible effect 
of abrupt variations of impedances between different 
media (e.g., between dura matter and cerebrospinal fluid). 
However, there is evidence that this may not be influential. 
First, our previous modeling work (Bedard et al., 2004) 
showed that abrupt variations of impedance have a 
negligible effect on low frequencies, suggesting that even 
in the presence of such abrupt variations should not play a 
role at low frequencies. Second, in the frequency range 
considered here, the skull and the skin are very close to be 
resistive at low frequencies (Gabriel et al., 1996b), so it is 
very unlikely that they play a role in the frequency scaling 
in EEG and MEG power spectra even at high frequencies. 

In the second part of the paper, we have analyzed 
simultaneous EEG and MEG signals recorded in four 
healthy human subjects while awake and eyes open (with 
desynchronized EEG). Because of the large number of 
channels involved, we used an automatic procedure 
(B-splines analysis) to calculate the frequency scaling. As 
found in previous studies (Pritchard, 1992; Freeman et al., 
2000; Bedard et al., 2006a), we confirm here that the EEG 
displays frequency scaling close to 1 // at low 
frequencie^]. However, this 1 // scaling was most typical 
of the midline channels, while temporal and frontal leads 
tended to scale with slightly larger exponents, up to 1/f 2 

5 Note that to compare scaling exponents between studies one 
must take into account that the electrode montage may influence the 
scaling. For example, in bipolar (differential) EEG recordings, if two 
leads are scaling as 1 / (A + /) and 1 / (B + /), the difference will have 
regions scaling as 1 /f 2 . 



(see Fig. H]\). The same pattern was observed in all four 
patients. 

This approach differs from previous studies in two 
aspects. First, in contrast to prior studies (such as Novikov 
et al., 1997; Linkenkaer-Hansen et al., 2001), we 
calculated the frequency scaling of all the sensors and did 
not confine our analysis to a specific region. Second, 
unlike other investigators (such as Hwa and Ferree, 
2002a,b), we did not limit our evaluations to either EEG 
or MEG alone, but rather analyzed the scaling of both 
type of signals simultaneously. Such a strategy enables us 
to provide an extended spatial analysis of the frequency 
scaling. It also provides a chance to compare the scaling 
properties of these signals in relation to their physical 
differences. 

For the MEG recordings, the frequency scaling at low 
frequencies was significantly lower compared to the EEG 
(see Fig. 0). This difference in frequency scaling was also 
accompanied by spatial variability patterns (see Fig. |4| 
showing three distinct regions: 1) a frontal area where the 
exponents had their highest values in the case of MEG; 2) 
a central area where the values of exponents of EEG and 
MEG get closer to each other and 3) a parietotemporal 
horseshoe region showing the lowest exponents for MEG 
with bimodal characteristics (Fig. [5]). In some cases, the 
scaling of the uncorrected and corrected MEG signal was 
also close to l/f, as reported previously (Novikov et al., 
1997). In the frontal area (FR mask), the scaling exponent 
of the EEG was generally larger. At Vertex (VX mask), 
EEG and MEG had similar values and at the 
Parietotemporal region (PT mask), MEG showed a 
bimodal property with a much broader range of scaling 
exponent in comparison to EEG (see Fig. |4|. Note that to 
avoid the effect of spurious peaks, Novikov et al. used the 
spectrum of signal differences and argued for the 
existence of a local similarity regime in brain activity. 
This approach fundamentally changes the spectral 
characteristics of Magnetometers (which measure the 
absolute magnitude of the magnetic induction) into a 
measure that only for the neighboring sensors 
approximates the behavior of the gradiometers (which 
measures the gradient of the magnetic induction). So it is 
not clear how to relate their values to the ones obtained 
here. 

To make sure that the differences of frequency scaling 
between EEG and MEG were not due to environmental or 
instrumental noise, we have used five different methods to 
remove the effect of noise. These methods are based on 
different assumptions about the nature and effect of the 
noise. A first possibility is to correct for the noise induced 
by the MEG sensors. It is known that the SQUID 
detectors used in MEG recordings are very sensitive to 
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environmental noise and they can produce 1 // noise 
(Hamalainen et al., 1993). Under this assumption, part of 
the scaling of the MEG could be due to "filtering" by the 
sensors themselves, which justifies a simple subtraction of 
scaling exponents to remove the effects of this filtering. 
Note that such empty-room recordings were not possible 
for the EEG, although the noise from the recording setup 
could be estimated (see Miller et al., 2009 for example). 
Because in some cases both MEG and emptyroom signals 
have similar scaling, a simple correction by subtracting 
the exponents would almost entirely abolish the frequency 
scaling while in other cases it may even revert the sign of 
the scaling exponent (see Fig.|4jl , Suppl. Fig. S3). 

However, if noise is not due to the sensors but is of 
additive uncorrelated nature, then another method for 
noise correction must be used. For this reason, we have 
used a second class of well-established methods 
consisting of spectral subtraction (Boll et al. 1979; Sim et 
al., 1998). Using three of such methods (LMSS, NMSS 
and WF) changed the scaling exponent, without 
fundamentally changing its spatial pattern (Fig. 3J)-F). 
The largest correction was obtained by non-linear 
methods which take into account the SNR information in 
the MEG signal. We also applied another class of method 
which uses the collective characteristics of all frequencies 
in noise correction (PLS). Similar to exponent subtraction, 
this method nearly abolished all the scaling of the MEG 
(Fig. Sfj). In conclusion, although different methods for 
noise subtraction give rise to different predictions about 
frequency scaling, all of the used methods enhanced the 
difference between EEG and MEG scaling. Thus, we 
conclude that the difference of EEG and MEG scaling 
cannot be attributed to noise, but is significant, therefore 
reinforcing the conclusion that the medium must be 
non-resistive. 

An alternative method to investigate this is the "Detrended 
Fluctuation Analysis" (DFA; see Kantelhardt et al., 2001; 
Linkenkaer-Hansen et al., 2001; Hwa and Ferree, 
2002a,b). Like many nonlinear approaches, DFA results 
are very vulnerable to the selection of certain parameters. 
Different filters severely affect the scaling properties of 
the electromagnetic signals to different extents, and 
therefore the parameters estimated through the DFA 
analysis could be false or lead to distorted interpretations 
of real phenomena (Valencia et al., 2008), and these 
effects are especially prominent for lower frequencies, 
which are precisely our focus of investigation here. One 
of the fields for which DFA can provide robust results is 
to analyze surrogate data with known characteristics. 
Although the use of DFA to evaluate the scaling 
exponents of EEG was vigorously criticized (Valencia et 
al., 2008), a previous analysis (Hwa and Ferree, 2002a,b) 
reported two different regions, a central and a more 



frontal, which somehow correlate with the FR and VX 
regions identified in our analysis. Similarly, a study by 
Buiatti et al. (2007) using DFA provided evidence for 
topographical differences in scaling exponents of EEG 
recordings. They report that scaling exponents were 
homogeneous over the posterior half of the scalp and 
became more pronounced toward the frontal areas. In 
contrast to Linkenkaer-Hansen et al. (2004) (where 
envelope of alpha oscillations was used for DFA 
estimation), this study uses the raw signal in its DFA 
analysis and yields values closer to those reported here. 

Both uncorrected signals and empty-room correction 
show that there is a fundamentally different frequency 
scaling between EEG and MEG signals, with near-1// 
scaling in EEG, while MEG shows a wider range at low 
frequencies. Although it is possible that non-neuronal 
sources affect the lower end (<lHz) of the evaluated 
frequency domain (Voipio et al., 2003), the solution to 
avoid these possible effects remain limited to invasive 
methods such as inserting the electrode into the scalp 
(Ferree et al., 2001) or using intracranial EEG recordings 
(similar to Miller et al., 2009). This approach would 
render wide range spatial recording as well as 
simultaneous invasive EEG and MEG recordings 
technically demanding or impractical. However, if 
technically feasible, such methods could provide a way to 
bypass non-neuronal effects at very low frequency. It 
could also provide a chance to evaluate the effects of 
spatial correlation on spectral structure at a multiscale 
level. 

The power spectral structure we observe here is consistent 
with a scenario proposed previously (Bedard et al., 
2006a): the 1 // structure of the EEG and LFP signals is 
essentially due to a frequency-filtering effect of the signal 
through extracellular space; this type of scaling can be 
explained by ionic diffusion and its associated Warburg 
impedancqj (see Bedard and Destexhe, 2009). It is also 
consistent with the matching of LFPs with multi-unit 
extracellular recordings, which can be reconciled only 
assuming a l/f filter (Bedard et al., 2006a). Finally, it is 
also consistent with the recent evidence from the transfer 
function calculated from intracellular and LFP recordings, 
which also showed that the extracellular medium is well 
described by a Warburg impedance (Bedard et al., 2010, 
submitted to this issue). If this non-resistive aspect of 
extracellular media is confirmed, it may influence the 
results of models of source localization, which may need 
to be reformulated by including more realistic 
extracellular impedances. 



6 Ionic diffusion can create an impedance known as the "Warburg 
impedance", which scales as 1/ i/fi), giving l/f scaling in the power 
spectra (Taylor and Gileadi, 1995; Diard et al., 1999). 
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In conclusion, the present theoretical and experimental 
analysis suggests that the scaling of EEG and MEG 
signals cannot be reconciled using a resistive extracellular 
medium. The 1 // structure of EEG with smaller scaling 
exponents for MEG is consistent with non-resistive 
extracellular impedances, such as capacitive media or 
diffusion (Warburg) impedances. Including such 
impedances in the formalism is non trivial because these 
impedances are strongly frequency dependent. The 
Poisson integral (the solution of Poisson's law 
V • D = —V • eW = p) would not apply anymore (see 
Bedard et al., 2004; Bedard and Destexhe, 2009). Work is 
under way to generalize the formalism and include 
frequency-dependent impedances. 

Finally, it is arguable that the scaling could also be 
influenced by the cancellation and the extent of spatial 
averaging of microscopic signals, which are different in 
EEG and MEG (for more details on cancellation see 
Ahlfors et al., 2010; for details on spatial sensitivity 
profile see Cuffin and Cohen, 1979). Such a possible role 
of the complex geometrical arrangement of underlying 
current sources should be investigated by 3D models 
which could test specific assumptions about the geometry 
of the current sources and dipoles, and their possible 
effect on frequency scaling. Such a scenario constitutes 
another possible extension of the present study. 



If we consider the differential expressions for the magnetic induction 
(Eq. l20t and electric potential (Eq. l2U . one can see that the frequency 
dependence of the ratio of their modulus is completely determined 
by the frequency dependence of the ratio of current density j'j! and 
j'f. In Appendix D, we show that this ratio is quasi-independent of 
frequency for a resistive medium, for low frequencies (smaller that 
~10 Hz), and if the current sources are of very low correlation. 

Thus, magnetic induction and electric potential can be very well 
approximated by: 

N 

Vf(r) = N < V >= N < £ SV l f > 

'=i (22) 
B/(r) = N < B >= N < £ SB' > 

for sufficiently small differential dendritic elements (N/l large). 

Because the functions of spatial and frequency are statistically 
independent, we can write the following expressions for the square 
modulus of the fields (see Eqs. l20l and l2U : 

\V f (r)\ 2 = N 2 \ < £ V'(r)G'/'(/) > | 2 = |V(r)| 2 |G(/)| 2 

l f (23) 
l|B/-(r)|| 2 = N 2 \\ < £ B'(r)G'/<(/) > || 2 = ||W(r)|| 2 |G(/)| 2 
i=i 

where G(/) =< Gf (/) >, V'(r) =< V'(r) > and W(r) =< 
B'(r) > . Thus, the scaling of the PSDs of the electric potential and 
magnetic induction must be the same for low frequencies (smaller 
than ~10 Hz) if the medium is resistive and when the current sources 
have very low correlation. 

B Differential expression for the magnetic induction 



Appendix 



A Frequency dependence of electric field and magnetic 
induction 



To compare the frequency dependence of magnetic induction and 
electric field, we evaluate them in a dendritic cable, expressed 
differentially. For a differential element of dendrite, in Fourier 
space, the current produced by a magnetic field (Ampere-Laplace 
law) is given by the following expression (see Appendix B): 



*B/(r) = £j?0O: 



Sv' 



(20) 



when the extracellular medium is resistive. Note that the source of 
magnetic induction is essentially given by the component of j'^ 
along the axial direction (J'r) within each differential element of 
dendrite because the perpendicular (membrane) current does not 
participate to producing the magnetic induction if we assume a 
cylindrical symmetry. 

For the electric potential, we have the following differential 
expression for a resistive medium (see Appendix C): 



1 SlUr 1 ) l j" f '(r') , 
/W ~4jry||r-i'|| ~ 4;ry ||r-r'|| 

where jj is the transmembrane current per unit of surface. 



(21) 



According to Maxwell equations, the magnetic induction is given by 

head 

where dv' — dx' l dx' 2 dx' 3 and 



dv' 



(24) 



||r-iT ||r-rT 
for a perfectly resistive medium. 

We now show that this expression is equivalent to Ampere-Laplace 
law. 

From the identity V x (gA) = g(V x A) + V'g x A, where 
V = e x + e y + e z -^, we can write: 



B ' (r) - S x <|F7j| > + ^ x V 'F^T ] dv ' (25) 

head 



Moreover, we also have the following identity 



(26) 



dhead 



where h is a unitary vector perpendicular to the integration surface 
and going outwards from that surface. Extending the volume integral 
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outside the head, the surface integral is certainly zero because the 
current is zero outside of the head. It follows that: 



head 

where dv' — dx l dx 2 dx 3 because 

1 , t 



r — r 



(27) 



r — r 



Eq.[^]is called the Ampere-Laplace law (see Eq. 13 in Hiimalainen 
et al., 1993). It is important to note that this expression for the 
magnetic induction is not valid when the medium is not resistive. 



We note that the expressions for the electric field and potential 
produced by each differential element of dendrite have the same 

frequency dependence because it is directly proportional to ^- for 
the two expressions. Also note that if the medium is resistive, then 
Yf = 7 and the frequency dependence of the electric field and 
potential are solely determined by that of the transmembrane current 

•m 
J f ■ 

D Frequency dependence of the ratio V/(x)/ jf (x). 



Finally, from the last expression, the magnetic induction for a 
differential element of dendrite can be written as: 



Sv' 



(28) 



For each differential element of dendrite, we consider the standard 
cable model, in which the impedance of the medium is usually 
neglected (it is usually considered negligible compared to the 
membrane impedance). In this case, we have: 



C Differential expression of the electric field and 
electric potential 



■ + ia>c m V, 



3V" 



,)Y<" 



(34) 



ff 



In this appendix, we derive the differential expression for the electric 
field. Starting from Eq. 1101 we obtain the solution for the electric 
potential: 



Vf(r) = ■ 



j_ m- v-j/(r') 



4^7/ JJJ ||r- 

head 



dv' 



(29) 



It follows that the electric field produced by the ensemble of sources 
can be expressed as: 



E/ (r) = -VV/(r) = 



1 



4iiYf 



V-j?(r': 



dV 



(30) 



such that every differential element of dendrite produces the 
following electric field: 



SE f (i) = — L -. m Sv' 

471:7/ ||r-r |P 



(31) 



The transmembrane current Slj^ obeys Slj^ — icopf(r')8v' because 
we are in a quasi-stationary regime in a differential dendritic element. 
Taking into account the differential law of charge conservation V • 
jy-(r') = — ifflpy(r'), we have: 



8E f (r) 



Slf(r') 



\TtYf 



4KY f 



,SS' 



(32) 



where /" is the density of transmembrane current per unit surface 
and 55' is the surface area of a differential dendritic element. This 
approximation is certainly valid for frequencies lower than 1000 Hz 
because the Maxwell-Wagner time (see Bedard et al., 2006b) of the 
cytoplasm (l2'r' = e/o ~ 10~ 10 s.) is much smaller than the typical 
membrane time constant of a neuron (t,„ ~ 5 — 20 ms). 

Finally the contribution of a differential element of dendrite to the 
electric potential at position r is given by 



«V,(r) = 



1 Slftf) _ 1 ]f(r>) 



AiiYf ||r-r'|| 47ty/ 1 1 r — r' 



SS' 



(33) 



where V", jK, c,„, r m et r, are respectively the membrane 
potential, the current density in the axial direction, the 
transmembrane current density, the specific capacitance (F/m 2 ), the 
specific membrane resistance (Q.nr) and the cytoplasm resistivity 
{Q.m). 



It follows that 



jf(x) r ; (l+;'(BT„,) dx 
where T„, = r m c m . 



[Zn(V/(x)] 



(35) 



Under in vivo-like conditions, the activity of neurons is intense and 
of very low correlation. This is the case for desynchronized-EEG 
states, such as awake eyes-open conditions, where the activity of 
neurons is characterized by very low levels of correlations. There is 
also evidence that in such conditions, neurons are in 
"high-conductance states" (Destexhe et al., 2003), in which the 
synaptic activity dominates the conductance of the membrane and 
primes over intrinsic currents. In such conditions, we can assume 
that the synaptic current sources are essentially uncorrelated and 
dominant, such that the deterministic link between current sources 
will be small and can be neglected (see Bedard et al., 2010). Further 
assuming that the electric properties of extracellular medium are 
homogeneous, then each differential element of dendrite can be 
considered as independent and the voltages V„, have similar power 
spectra. 



In such conditions, we have: 

y;"(x) = F'"(x)G" , (/) 



(36) 



Note that this expression implies that we have in general for each 
differential element of dendrite: 



f jf(x) = F'"(x)(i±if-)G»'(/) 

1 //'(*) = 4^G"'(/) = F''(x)G'"(/) 
according to Eq. 1341 



(37) 
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It follows that 

4t4 = -r^ 1 7 ~[ln(F(x))] » — • 4-[/n(F(x))l (38) 

Thus, for frequencies smaller than l/(o)T m ) (about 10 to 30 Hz for 
T m of 5-20 ms), the ratio .{„ ^ will be frequency independent, and for 
each differential element of dendrite, we have: 

r/;.-(x) = f"'(x)g'"(/) 

\ ;>(x) = F'^G'n/) ^' 
for frequencies smaller than ~10 Hz. 
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Supplementary material 



spectrum (original signal) and noise magnitude spectrum 
estimate ("noise"), respectively, k is the frequency index, 
while a* and bk are linear coefficient parameters of the 
summation. Spectral subtraction methods fall into three 
main categories (Sim et al., 1998). The simplest of all, a 
linear method where = bk = 1, 0£=2, following Boll et 
al. (1979) was used here. This linear multiband spectral 
subtraction (LMSS) method is well-established for noise 
subtraction (see Loizou, 2007 for a comparative study of 
noise subtraction methods). 



Supplementary methods 



We give details below to some of the methods and 
quantities used in the Results. 



SNR 



Two of the used methods for noise-correction are based 
on band-specific signal-to-noise ratio (SNR) in order to 
cancel the effects of background colored-noise in the 
spectra of interest. In each subject, average PSD was used 
to calculate signal-to-noise ratio (SNR). For SNR 
calculation, few frequency bands were defined based on 
the categorization in Buzsaki & Draguhn (2004): 0-10 Hz 
(Slow, Delta and Theta), 11-30 Hz (Beta), 30-80 Hz 
(Gamma), 80-200 Hz (Fast oscillation), 200-500 Hz 
(Ultra-fast oscillation). SNR was calculated as: 



SNR bi 



(40) 



for a given band "b" and sensor "i", "n" is the frequency 
resolution of that band. This method was applied on 
individual average PSD as well as shape preserving spline 
of each average PSD where each PSD was fist smoothed 
in loglO scale using a shape preserving spline, i.e, 
Piecewise Cubic Hermite Interpolating Polynomial 
(PCHIP). 



An improved method, with = 1 and bk = v, where "v" is 
the oversubtraction factor. This method uses 
oversubtraction and introduces a spectral flooring to 
minimize residual noise and musical noise (Berouti et al., 
1979). A second category of spectral subtraction is based 
on elk = bk =f(k). Third and the most robust methods are 
based on a non-linear multiband subtraction (NMSS) 
where cif, = 1 and bk = v(k); i.e., the oversubtraction factor 
is adjusted based on a specfic band's SNR. These methods 
proposed by (Kamath and Loizou, 2002; Loizou, 2007) 
are suitable for dealing with colored noise (Boubakir et 
al., 2007; Sim et al., 1998), a case similar to MEG 
recordings. The spectrum is divided into N 
non-overlapping bands, and spectral subtraction is 
performed independently in each band. The Eqs. [41] is 
simply reduced to: 



\S^k)\ 2 = \Yi{k)\ 2 - a,«S, \D^k)\ z ,bi <k<, 



(42) 



where bj and e, are the beginning and ending frequency 
bins of the ith frequency band, a,- is the overall 
oversubtraction factor of the ith band and 5,- is a tweaking 
factor. The band specific oversubtraction factor a,- is a 
function of the segmental SNR, of the ith frequency band. 
After calculating bandspecific SNR (Eqs.l40l, we used the 
product of lower 10 percent of crosssubject average SNR 
and standard deviation of SNRj to estimate the a, (5, 
subtraction coefficient. Next, simply by multiplying the 
noise PSD by this coefficient and subtracting it from the 
measured PSD, the enhanced PSD was achieved. 



Multiband spectral subtraction 

Assuming the additive noise to be stationary and 
uncorrected with the clean signal, nearly most spectral 
subtraction methods can be formulated using a parametric 
equation: 

\sJk)\ a =a k \Y(k)\ a -b k \D(k)\ a (41) 

where |S*|, \Yk\ and |Z)#| refer to enhanced magnitude 
spectrum estimate (corrected signal), the noisy magnitude 



Wiener filter (WF) spectral enhancement 

The principle of the Wiener filter is to obtain an estimate 
of the clean signal from that of the noisy measurement 
through minimizing the Mean Square Error (MSE) 
between the desired and the measured signal (Lim et al., 
1979; Abd El-Fattah et al., 2008). In the frequency 
domain, this relation is formulated as filtering transfer 
function: 

WF ^=mrm (43) 
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where, as before, P s (k) and P„(k) refer to enhanced power 
spectrum estimate and noise power spectrum estimate 
respectively for a signal frame and k is the frequency 
index. Based on the definition of SNR as, the ratio of 
these two elements, one can formulate the WF as: 



A. Mean and standard deviation 



WF K = [1 



1 



SNR k ' 



(44) 



After calculation of bandspecific WF, the noisy signal is 
simply muliplied by the WF to obtain the enhanced signal. 



Partial least square (PLS) approximation of non-noisy 
spectrum 



Partial least squares (PLS) regression, combines 
"Principal component analysis" (PCA) and "Multiple 
linear regression" (Abdi, 2010; Abdi and Williams, 2010). 
While PCA finds hyperplanes of maximum variance 
between the response and independent variables, PLS 
projects the predicted variables and the observable 
variables to a new space. Then from this new space, it 
finds a linear regression model for the projected data. 
Next, using this model, PLS finds the multidimensional 
direction in the X space that explains the maximum 
multidimensional variance direction in the Y space (Abdi, 
2010; Garthwaite, 1994). If X is the PSD of noise 
measurement and Y is the PSD of the measured signal 
contaminated with background noise, one can use PLS to 
"clean" one matrix (Y) by predicting Y from X and then 
using the residual of the prediction of Y by X as the 
estimate of pure PSD. The patterns of the awake spectrumS 
that statistically resembles the patterns of emptyroomj? 
spectral noise are those that should be removed. As during 
PLS algorithm, the data is mean subtracted and 
z-normalized, the predection of Y from X is an 
approximate of the zscored PSD. Therefore, the reseidual 
Y, which is taken as the spectral features that can not be 
predicted by noise, also has zscored values. It has too be 
emphasized that this approach of denoising only works in 
the spectral but not the time domain. 
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Supplementary figures 





EEG 


MEG (awake) MEG(empty) LMSS 


All 


-1.33 ± 0.19 


-1.24 ± 0.26 


-1.04 ±0.13 -1.24 ±0.28 


FRROI 


-1.36 ±0.25 


-0.97 ± 0.10 


-0.97 ± 0.06 -0.96 ±0.11 


VXROI 


-1.21 ±0.13 


-1.36 ±0.10 


-1.10 ± 0.09 -1.36 ±0.10 


PTROI 


-1.36 ±0.12 


-1.30 ±0.29 


-1.08 ±0.15 -1.31 ±0.32 




NMSS 


WF 


PLS 


ES 


All 


-1.06 ±0.29 


-1.05 ±0.27 


-0 SO + 1 1 


-0 20 + 23 


FRROI 


-0.76 ± 0.09 


-0.76 ± 0.08 


-0 40 + 05 


-0 00 + 09 

\J . \J\J _L_ U.U7 


VXROI 


-1.14 ±0.11 


-1.12 ±0.11 


-0.50 ± 0.04 


-0.26 ± 0.08 


PTROI 


-1.16 ±0.32 


-1.14 ±0.30 


-0.54 ±0.11 


-0.22 ± 0.26 


B. Pearson correlation of EEG vs. 








MEG LMSS 


NMSS 


WF PLS 


ES 


All 


0.29 0.29 


0.32 


0.33 0.37 


0.35 


FRROI 


0.41 0.39 


0.32 


0.37 0.01 


0.17 


VXROI 


-0.17 -0.10 


-0.15 


-0.13 0.01 


-0.28 


PTROI 


0.35 0.34 


0.38 


0.39 0.46 


0.41 


C. Kendall Rank Corr of EEG vs. 








MEG LMSS 


NMSS 


WF PLS 


ES 


All 


0.21 0.21 


0.24 


0.25 0.29 


0.23 


FRROI 


0.29 0.23 


0.21 


0.27 -0.06 


0.12 


VXROI 


-0.03 0.04 


-0.04 


-0.03 0.07 


-0.09 


PTROI 


0.23 0.23 


0.26 


0.26 0.30 


0.27 



Table 2 ROI statistical comparison for different noise correction 
methods. A. mean and std of frequency scale exponent for all regions 
and individual ROI. B. numerical values of linear Pearson correlation. 
C. rank-based Kendall correlation. 



MAG 



GRAD1 



GRAD2 






Figure Frequency spectra of magnetometers and gradiometers. 
Comparison of awake (blue) vs empty-room (red) recordings 
between Magnetometers (MAG) and Gradiometers (GRAD1, 
GRAD2) in a sample subject. As for the EEG, the MEG signal is 
characterized by a peak at around 10 Hz, which is presumably due 
to residual alpha rhythm (although the subject had eyes open). This 
is also visible from the MEG signals (Fig. [JJ as well as from their 
PSD (Fig. [5] and MAG panel here). The power spectrum from the 
empty-room signals also show a peak at around 10 Hz, but this 
peak disappears from the gradiometer empty-room signals, while the 
10 Hz peak of MEG still persists for gradiometers awake recordings. 
This suggests that these two 10 Hz peaks are different oscillation 
phenomena. All other subjects showed a similar pattern. 
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Figure S2: Signal-to-noise ratio (SNR) of Magnetometers (MAG) for 
multiple frequency bands: 0-10 Hz (Slow, Delta and Theta), 11-30 
Hz (Beta), 30-80 Hz (Gamma), 80-200 Hz (Fast oscillation), 200- 
500 Hz (Ultra-fast oscillation). In the scatterplots, red astrisks relate 
to individual sensors and the blue line is the band-specific mean 
across the sensors. In boxplots, the box has lines at the lower quartile, 
median (red), and upper quartile values. Smallest and biggest non- 
outlier observations (1.5 times the interquartile range IRQ) are shown 
as whiskers. Outliers are data with values beyond the ends of the 
whiskers and are displayed with a red + sign. In all subjects, the 
SNR shows a band-specific trend and has the highest value for lower 
frequencies and gradually drops down as band frequency goes up. As 
the frequency drops, the variability of SNR (among sensors) rises; 
therefore, the SNR of the lowest band (1-10 Hz) shows the highest 
sensors-to-sensor variability and the highest SNR in comparison to 
other frequency bands. 
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Awake LMSS NMSS PLS ES 

Figure S3: Noise correction comparison. Every horizontal line 
showes a voxel of the topographical maps shown in Fig. [4] sorted 
based on the scaling exponent values of awake MEG (left stripe). 
Using a continuous color spectrum, these stripes show that minimal 
correction is achived by LMSS. As indicated in the text, the 
performance of this method is not reliable due to the nonlinear 
nature of SNR (see Suppl. Fig. S2). NMSS yields higher degree of 
correction. WF performs almost identical to NMSS (not shown here). 
Exponent subtraction almost abolishes the sacling all together (far 
right stripe). PLS results in values between NMSS and "Exponent 
subtraction". For details of each of these correction procedures, see 
Methods. LMSS, NMSS and WF rely on additive uncorrected nature 
of noise. "Exponent subtraction" assumes that the noise is intrinsic to 
SQUID. PLS ascertains the characteristics of noise to the collective 
obeserved pattern of spectral domain across all frequencies. See text 
for more details. 



